load("~/Documents/MiGASti/Databases/gene_matrix.RData")
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
#exclude samples that did not pass QC filtering
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
genes_counts_filt = genes_counts[, !colnames(genes_counts) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_counts_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)
gencode_30 = read.table("~/Documents/MiGASti/Databases/ens.geneid.gencode.v30")
colnames(gencode_30) = c("ensembl","symbol")
#order metadata and genes counts
genes_counts_ordered <- genes_counts_filt[,rownames(metadata_filt)]
#head(genes_counts_ordered)
all(rownames(metadata_filt) == colnames (genes_counts_ordered)) #TRUE
# Counts, TPM, voom and metadata filtered 
# 496 samples!  
# Filter for genes not expressed: 50 % of the samples! 
# dim(genes_counts_ordered) #18997   496
# dim(metadata_filt) #496  38

Preparing the samples for DEG

#remove uncultured samples
metadata_cultured <- metadata_filt[metadata_filt$Stimulation != "ununstim",]
#check numbers
#dim(metadata_cultured)
#check numbers per stimulation
#table(metadata_filt$Stimulation)
#select only THA samples in metadata
metadata_THA = metadata_cultured[metadata_cultured$Region=="THA",]
#select only THA samples in genes counts
genes_counts_THA <- genes_counts_ordered[,metadata_THA$Sample]
#order metadata and genes_counts
genes_counts_THA_ordered <- genes_counts_THA[,rownames(metadata_THA)]
#check ordering
all(rownames(metadata_THA) == colnames (genes_counts_THA_ordered)) #TRUE

[1] TRUE

#round counts; deseq2 can only handle integers
genes_counts_THA_ordered <- round(genes_counts_THA_ordered, digits=0)

#make sure covariate variables are the right format 
#cannot create dds object with numeric values
metadata_THA$Donor_id <- as.factor(metadata_THA$Donor_id)
metadata_THA$age <- as.integer(metadata_THA$age)
metadata_THA$sex <- as.factor(metadata_THA$sex)
metadata_THA$Stimulation <- as.factor(metadata_THA$Stimulation)
metadata_THA$picard_pct_ribosomal_bases = scale(metadata_THA$picard_pct_ribosomal_bases)
metadata_THA$picard_pct_mrna_bases = scale(metadata_THA$picard_pct_mrna_bases)
metadata_THA$picard_pct_pf_reads_aligned = scale(metadata_THA$picard_pct_pf_reads_aligned)
metadata_THA$picard_pct_percent_duplication = scale(metadata_THA$picard_pct_percent_duplication)

#adjust for: ~ age + (1|donor_id) + picard_pct_ribosomal_bases + picard_pct_mrna_bases +   picard_pct_percent_duplication + picard_pct_pf_reads_aligned

DESeq2 of THA samples

#createDeSEQ2 object for LPS
dds <- DESeqDataSetFromMatrix(countData = genes_counts_THA_ordered,
                              colData = metadata_THA,
                              design = ~ age + sex + picard_pct_ribosomal_bases + picard_pct_mrna_bases + picard_pct_percent_duplication + picard_pct_pf_reads_aligned + Stimulation) 
#variable of interest at end of the formula

#Make sure that control group is set as the reference group
dds$Stimulation <- relevel(dds$Stimulation, ref="unstim")
table(dds$Stimulation)

unstim IFNy LPS 22 9 22

#head(dds)

#filter: CPM > 1 in 50% of the samples 
keep.exp = rowSums(cpm(genes_counts_THA_ordered) > 1) >= 0.5*ncol(genes_counts_THA_ordered)
dds = dds[keep.exp,]

#Run differential expression 
dds <- DESeq(dds, betaPrior = FALSE)
resultsNames(dds)

[1] “Intercept” “age”
[3] “sex_m_vs_f” “picard_pct_ribosomal_bases”
[5] “picard_pct_mrna_bases” “picard_pct_percent_duplication” [7] “picard_pct_pf_reads_aligned” “Stimulation_IFNy_vs_unstim”
[9] “Stimulation_LPS_vs_unstim”

DESeq2: LPS vs unstim

Number of differentially expressed genes

# generate results table for LPS vs unstim
res_LPS <- results(dds, name="Stimulation_LPS_vs_unstim")
sum(res_LPS$padj < 0.05, na.rm=TRUE)

[1] 190

resOrdered_LPS <- res_LPS[order(res_LPS$pvalue),] 
resOrdered_LPS <- as.data.frame(resOrdered_LPS)

Volcano plot LPS vs unstim

head(res_LPS)

log2 fold change (MLE): Stimulation LPS vs unstim Wald test p-value: Stimulation LPS vs unstim DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue ENSG00000000003.14 17.9465 0.4108771 0.277869 1.478672 0.139228 ENSG00000000419.12 297.9371 -0.2043851 0.232001 -0.880966 0.378336 ENSG00000000457.14 101.2246 0.1552885 0.169181 0.917886 0.358679 ENSG00000000460.17 37.0403 -0.0850284 0.234128 -0.363170 0.716478 ENSG00000000938.13 574.9011 -0.1063089 0.212124 -0.501163 0.616256 ENSG00000000971.15 41.6091 0.7399651 0.353692 2.092117 0.036428 padj ENSG00000000003.14 NA ENSG00000000419.12 0.662735 ENSG00000000457.14 0.648166 ENSG00000000460.17 0.878170 ENSG00000000938.13 0.822562 ENSG00000000971.15 0.277225

with(res_LPS, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2)))
with(subset(res_LPS, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))

MA plot LPS vs unstim

#The function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.

plotMA(res_LPS, ylim=c(-2,2))

### TOP differentially expressed genes LPS vs unstim

setDT(resOrdered_LPS, keep.rownames = "ensembl")
resOrdered_LPS <- left_join(resOrdered_LPS, gencode_30, by = "ensembl")
resOrdered_LPS_top = resOrdered_LPS[order(resOrdered_LPS$padj) ,]
setDT(resOrdered_LPS_top, keep.rownames = "ensembl")
resOrdered_LPS_top = resOrdered_LPS_top[, c("ensembl", "symbol", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")]
createDT(resOrdered_LPS_top)
write.table(resOrdered_LPS_top, "DEG_LPS_THA.txt")

DESeq2: IFNy vs unstim

Number of differentially expressed genes

# generate results table for IFNy vs unstim
res_IFNy <- results(dds, name="Stimulation_IFNy_vs_unstim")
sum(res_IFNy$padj < 0.05, na.rm=TRUE)

[1] 43

resOrdered_IFNy <- res_IFNy[order(res_IFNy$pvalue),] 
resOrdered_IFNy <- as.data.frame(resOrdered_IFNy)

Volcano plot IFNy vs unstim

head(res_IFNy)

log2 fold change (MLE): Stimulation IFNy vs unstim Wald test p-value: Stimulation IFNy vs unstim DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue ENSG00000000003.14 17.9465 0.362369 0.350506 1.033846 0.301208 ENSG00000000419.12 297.9371 -0.263876 0.303501 -0.869438 0.384607 ENSG00000000457.14 101.2246 0.304341 0.219580 1.386015 0.165742 ENSG00000000460.17 37.0403 -0.211603 0.304589 -0.694716 0.487233 ENSG00000000938.13 574.9011 -0.175792 0.277655 -0.633130 0.526649 ENSG00000000971.15 41.6091 0.617750 0.461793 1.337719 0.180988 padj ENSG00000000003.14 NA ENSG00000000419.12 0.999854 ENSG00000000457.14 0.999854 ENSG00000000460.17 0.999854 ENSG00000000938.13 0.999854 ENSG00000000971.15 0.999854

with(res_IFNy, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-10,10)))
with(subset(res_IFNy, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))

MA plot IFNy vs unstim

plotMA(res_IFNy, ylim=c(-2,2))

TOP differentially expressed genes IFNy vs unstim

setDT(resOrdered_IFNy, keep.rownames = "ensembl")
resOrdered_IFNy <- merge(resOrdered_IFNy, gencode_30, by = "ensembl")
resOrdered_IFNy_top = resOrdered_IFNy[order(resOrdered_IFNy$padj) ,]
setDT(resOrdered_IFNy_top, keep.rownames = "ensembl")
resOrdered_IFNy_top = resOrdered_IFNy_top[, c("ensembl", "symbol", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")]
createDT(resOrdered_IFNy_top)
write.table(resOrdered_IFNy_top, "DEG_IFNy_THA.txt")

Create genelists logFC > 1 or <-1

resOrdered_LPS_p <- subset(resOrdered_LPS_top, padj < 0.05)
resOrdered_LPS_LFC <- subset(resOrdered_LPS_p, log2FoldChange > 1 | log2FoldChange < -1)
write.table(resOrdered_LPS_LFC, "LPS_THA_FC1.txt")

resOrdered_IFNy_p <- subset(resOrdered_IFNy_top, padj < 0.05)
resOrdered_IFNy_LFC <- subset(resOrdered_IFNy_p, log2FoldChange > 1 | log2FoldChange < -1)
write.table(resOrdered_IFNy_LFC, "IFNy_THA_FC1.txt")